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; ABSTRACT 

^ f We explore the dynamical signatures imprinted by baryons on dark matter haloes 

during the formation process using the Overwhelmingly Large Simulations (OWLS), 
a set of state-of-the-art high-resolution cosmological hydrodynamical simulations. We 
present a detailed study of the effects of the implemented feedback prescriptions on 
I the orbits of dark matter particles, stellar particles and subhaloes, analysing runs 

with no feedback, with stellar feedback and with feedback from supermassive black 
holes. We focus on the central regions (0.25r2oo) of haloes with virial masses 6 x 
^! 10^^{r^ 7 X lO")/i"iM0 a.t z = 0(2). We also investigate how the orbital content 

(relative fractions of the different orbital types) of these haloes depends on several key 
parameters such as their mass, redshift and dynamical state. The results of spectral 
■ analyses of the orbital content of these simulations are compared, and the change 

Ti^lj- ■ in fraction of box, tube and irregular orbits is quantified. Box orbits are found to 

dominate the orbital structure of dark matter haloes in cosmological simulations. 
There is a strong anticorrelation between the fraction of box orbits and the central 
baryon fraction. While radiative cooling acts to reduce the fraction of box orbits, strong 
feedback implementations result in a similar orbital distribution to that of the dark 
matter only case. The orbital content described by the stellar particles is found to be 
remarkably similar to that drawn from the orbits of dark matter particles, suggesting 
that either they have forgotten their dynamical history, or that subhaloes bringing in 
j^j ■ stars are not biased significantly with respect to the main distribution. The orbital 

I content of the subhaloes is in broad agreement with that seen in the outer regions of 

the particle distributions. 

Key words: galaxies: haloes - galaxies: clusters: general - galaxies: evolution - galax- 
ies: kinematics and dynamics - methods: numerical - cosmology: theory 
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1 INTRODUCTION 

Dark matter structure formation is well understood within 
the standard cosmological model. Haloes are thought to form 
hierarchically, through the merging and accretion of smaller 
systems. As such, there should be observational signatures of 
these merging processes in the resulting remnants, providing 
dynamical information about their formation histories. We 
investigate the orbital content of dark matter haloes in order 
to explore what signatures may result. 
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Dark matter haloes formed in a Acold dark mat- 
ter (ACDM) cosmology appear to share a nearly uni- 
versal internal morphology; they have density profiles 
that are well des cribed by the Navarro, Frenk & White 
(hereafter NFW; Il996l . 1 19971 ') profile and pseudo-phase 
space densities with a constant (a ~ 1.9) power-law 
slope llTavlor fc Navarroll200ll:lDehnen fc McLaughlinll2005l : 



iBarnes et all l2006l : iLudlow et al.l boill ). There 



IS a uni- 



versal relation between the radial density profile slope 
and the velocity a nisotropy within the i n ner region of 
dark matter haloes (|Hansen fc Moor3 l2006l : iNavarro et al.l 
I2OI0I ) and the velocity distribution function is found to 
have a universal shape ([Hansen et al.l |2006| V Dark mat- 
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ter haloes are thought to have spin distributions that 
are reasonably well characterised by a log-normal dis- 



tribution ( 


Bullock et al.' '200l': 'Bailin & Steinmetd l2005l: 


Bett et all 


2007,: .Maccio et al.. ,200^) and are thouE 


fht to 


be triaxial (IFrenk et alj 19881: Dubinski & Carlberd 


1991; 


Warren et all 19921: Cole & Lacev 19961: Jins 


& Sutd 


2002: 


Baihn & Steinmetz 20051: Allgood et al.l 2006 


: Maccio et al. 


20061: Bett et al. 20071: Jeeson-Daniel et all 2011^. Here we 



explore the underlying orbital distribution of dark matter 
haloes evolved within the standard cosmological model. 

Of course, a complete understanding of structure for- 
mation requires careful consideration of the baryonic com- 
ponent, particularly in the inner region of haloes. The 
baryonic physics processes involved are largely uncertain 
and have become the focus of galaxy formation stud- 
ies. It is, however, well established that the condensa- 
tion of baryons to the centres of dark matter haloes 
tends to result in the halo becoming more spherical or 
axisymmetric, and that this is a direct consequence of 

the impact of baryons on the orbital content of the 
, , , „ 1 ^ — Z — T~o — TT^ 1 1. „otr' I.- — T> 
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1996; 
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Jesseit et al. 20051: IDebattista et all I2OO8I: IValluri et al.l 



varying implementations of baryonic physics to systemati- 
cally explore the effect of baryons on the orbital content of 
the haloes. We do so by performing a spectral analysis on 
the snapshots of the simulation, assuming that the potential 
remains static during an orbit. 

Orbits are divided into groups, or families, according 
to the phase space they cover. Families of regular orbits 
have similar morphologies because they conserve similar in- 
tegrals of motion. These isolating integrals restrict the re- 
gion of phase space available to an orbit. Each lowers, by 
one, the dimensionality of the region available to the or- 
bit. As such, an orbit is shaped by its isolating integrals. 
Axisymmetric potentials have two classical integrals of mo- 
tion: energy E and the ^-component of the angula r momen- 
tum Lz a nd a third non-classical isolating integral (|Lindbladl 
I1933I . [Con topoulos 1960, Binncy & Spcrgcl 198^). However, 
it is wel l-known that man y elliptical galaxies are not axisym- 
metric l|Franx et al.ll 19911) . In more relevant triaxial systems, 
there is a second non-classical integral that i s likely to play 
a rol e in shaping the structure of the system (jSchwarzschildl 
11973 1. For this reason numerical simulations provide an im- 
portant tool in understanding the orbital content of such 
systems. 

Within simple generic triaxial models, regular orbits 
can be divided into two main families; box and tube orbits. 
Tube orbits are further divided accord ing to their orienta- 
tion into ma jor- and minor-axis tubes (|Schwarzschildlll979l : 
IStatled[l987l ). Box orbits are free to pass close to the centre 
of the potential and their orbit-averaged angular momentum 
is zero. While box orbits show no sense of rotation, tube or- 
bits tend to rotate around the centre of the system, avoiding 
the centre. Box orbits are fundamentally important, as they 
are thought to be responsible for conveying information from 
the central regions of a halo to the outer parts of the system 
and are required to support the triaxial halo. 

Several authors have investigated the orbital content of 
analytic potentials and remnants of simulated disc mergers. 
An initial attempt to classify the orbital content of simu- 



lated merger remnants was conducted bv lBarned l| 19921 ) who 
simulated a small sample of merging encounters between 
equal mass disc galaxies. Orbits were classified according to 
changes in the sign of the angular momentum vector. The 
shapes and kinematic properties of the remnants were found 
to be related to the initial spin vectors and other encounter 
paramet ers. By including gas dyna mics in these merging 
galaxies, iBarnes fc HernquistI (|l996h showed the dramatic 
effect gas can have on the structure of the resulting rem- 
nant. Torques experienced during the merger act to remove 
angular momentum from the gas, causing it to flow inwards 
to form a central mass concentration. They found the depth 
of the potential well to be highly correlated with the stellar 
kinematics and that gas ac ts to destabilise box orbits (as 
discussed in lDubinskilll993 ). This causes minor-axis tubes 
to become dominant and results in a more oblate remnant. 

I Jesseit et al.l l|2005h studied a statistical sample of disc 
ga laxy mergers, using th e auto mated spectral classification 
of ICarpintero fc AguilaJ (|l998t ) to quantify the orbital con- 
tent of the resulting remnants. They found that the most 
abundant orbital classes were box and minor-axis tube or- 
bits. While the inner regions of the simulated remnants were 
dominated by box orbits, tube orbits became more impor- 
tant at intermediate radii. They also found that the ratio of 
these two classes of orbits played a role in determining the 
basic properties of the remnant. Minor-axis-tube-dominated 
haloes were found to be discy, while those dominated by box 
orbits were boxy. Major-axis tubes were found to be domi- 
nant in prolate remnants. Again, it was noted that gas af- 
fects the fraction of box orbits, causing an increase in the 
population of minor-axis tubes. 

IDebattista et all l|2008l ) studied the impact of growing 
a central disc on the orbital content of a halo. They find 
that while the central concentration does result in rounder, 
more radially anisotropic haloes, the halo's shape is essen- 
tially returned to its original state if the disc is artificially 
'evaporated'. This indicates that the character of the or- 
bits is not generally changed by the central mass concentra- 
tion, the box orbits are not destroyed but simply become 
ro under in line with the potential. This is also considered 
m IVaUuri et all (|2010D who explore the orbital evolution in- 
duced by baryonic condensation in triaxial haloes. They find 
that the evolution depends on the radial distribution of the 
baryonic component, and that a massive compact central 
mass will result in the scattering of a large fraction of both 
box and long-axis tube orbits even at fairly large pericentric 
distances. 

A comprehensive study of the orbital structur e of 1: 1 
merger remnants can be found in iHoffman etahl HoiO). 
Mergers between equal mass discs at varying initial gas frac- 
tions (ranging from to 40 per cent) were simulated, tak- 
ing into account both star formation and feedback. They 
showed that, by varying the fraction of gas in a merger, a 
wide range of kinematic structures can be produced. The 
remnants formed in these simulations are typically prolate- 
triaxial. The central regions are dominated by box orbits, 
while tube orbits dominate further out. The inclusion of gas 
acts to decrease the fraction of stellar particles on box or- 
bits in the central region, replacing them with minor-axis 
tubes. The remnants were found to become progressively 
more oblate as the gas fraction is increased. Outside of 1.5 
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Re (where Re is the effective radius) the remnants are found 
to be largely unaffected by the addition of gas. 

This work aims to extend the previous work on the effect 
of baryons on orbital structure by comparing several models 
for the feedback implementation within realistic cosmologi- 
cal simulations. 

The outline of the paper is as follows. In Section [2] we 
briefly review the simulations used for this study, describ- 
ing the different baryonic physics implemented and the halo 
sample extracted for analysis. In Section [3] we describe the 
method used to define the orbital content of the haloes. Our 
main results are presented in Section [4] and we finish with 
a summary of our conclusions in Section [S] Resolution is- 
sues/convergence tests and the effect of halo definition are 
discussed in Appendix A. 



2.1 Baryon physics 

In this work we have analysed a subset of five of the OWLS 
runs to explore the effect of baryons on the orbits of dark 
matter haloes by studying simulations with varying levels of 
feedback. We have considered a dark matter only run; a run 
which includes baryons and primordial element line cooling 
but no feedback; a weak stellar feedback run; a strong stellar 
feedback run and a run which includes feedback from super- 
massive black holes in addition to weak stellar feedback. The 
simulations used and their implemented subgrid physics are 
summarized in Table [T] The baryon runs consider gas cool- 
ing and star formation as well as feedback from stars and 
active galactic nuclei (AGN). The implementation of each 
of these processes is discussed briefly below. 



2 NUMERICAL SIMULATIONS 

The haloes used for this analysis were extracted from the 
Overwhelmingly Large Simulations (OWLS). OWLS is a set 
of high-resolution cosmological simulations run with vary- 
ing implementations of the subgrid physics. For detailed in- 
formation about th ese simulations we refer the reader to 
ISchave et al] (|201(]l l. Here we discuss only briefly the perti- 
nent details of the subset of simulations used for our analy- 
sis. The cosmological parameters were taken from the third 
year Wilkinson Microwav e Anisotropy Probe (WMAP3) re- 
suhs l|Spergel et al.ll2007l ). with: fi^ = 0.238, = 0.762, 
Q.b = 0.0418, h = 0.73, n = 0.95 and as= 0.74. The pri- 
mordial baryonic mass fraction of hydrogen (helium) is as- 
sumed to be 0.752 (0.248). For all runs, cosmological ini- 
tial conditions were constructed using a linear power spec- 
trum based on a tr ansfer function generated with CMBFAST 
l|Seliak fc ZaldarrTaga 1996 ). The initial posit ions and ve - 
locities w ere computed from a glass-like state (|Whitelll996l ) 
usmg the lZel'dovichl |l973) approximation. 

The sim ulations were r un using a modified version of 
GADGET-3 (ISpringell l2005l : ISpringel et all |2008| ) to foUow 
the evolution of 512^ dark matter particles and 512^ gas 
particles in cubes of comoving lengths 25 and 100/i~^Mpc. 
In the dark matter only run the particle mass is 7.65 x 10^ 
and 4.93 x 10* /i'^Mq in the 25 and 100 /i"^ Mpc boxes 
respectively. The baryon runs follow the baryonic compo- 
nent with smooth particle hydrodynamics (SPH), where the 
number of neighbours A'^ngb for the SPLI interpolation was 
set to 48. The mass of the particles in the baryon runs is di- 
vided between the gas and dark matter particles according 
to the universal baryon fraction, J^"'" = f^t/ilm = 0.176, 
such that the dark matter (gas) mass in the 100 /i~^Mpc run 
is 4.06(0.87) X 10^ h'^MQ and 6.34(1.35) x 10*^ /i-^M© for 
the 25 /i~^Mpc box. Baryonic particle masses are allowed to 
change during the simulation due to mass transfer from star 
to gas particles. The comoving gravitational force softening 
was set to 1/25 of the initial mean interparticle spacing but 
is limited to a maximum physical scale of 2 (0.5) ^~^kpc for 
the 100 (25) /i~^Mpc boxes. For a given box size the same 
initial conditions are used in each run. This allows us to fol- 
low the same haloes with different implementations of the 
subgrid physics and provides a unique opportunity to sys- 
tematically study the effects of the subgrid physics on haloes 
evolved within a cosmological framework. 



2.1.1 Cooling 

The cooling rates are computed element-by-element in 
the presence of the co smic microwave background and 
the llfaardt fc Madaul (|200ll ') model for the evolving 
ultraviolet/X-ray background radiation from quasars and 
galaxies. Contributions from hydrogen, helium, carbon, ni- 
trogen, oxygen, neon, magnesium, silicon, sulphur, calcium 
and iron are considered. The contributions are interpolated 
as a function of density, te mperature and r edshift from pre- 
computed CLOUDY tables IjWiersma et al ]|2009), assummg 
the gas to be optically thin and in photo-ionisation equilib- 
rium. In the NOSN_NOZCOOL run, cooling rates are calcu- 
lated using primordial element abundances. Cooling by both 
Bremsstrahlung emission and Compton cooling via interac- 
tions between the gas and cosmic microwave background is 
also taken i nto account. Reion i sation is modelled by 'switch- 
ing on' the iHaardt fc Madaul (|200ll ) background at z = 9. 
Collisional equilibrium is assumed before reionisation, and 
photo-ionisation after z = 9. 



2.1.2 Star formation and evolution 

Star formation is modelled by converting gas particles 
into collisionless stellar particles (representing a sim- 
ple stellar population) according to the prescription of 
ISchave fc Dalla Vecchial (|2008l ). A star formation density 
threshold of uh > 0.1 cm~^ is adopted; above this den- 
sity an effective equation of state (P cx p'^'^^) is imposed, 
where 7ofi is set to 4/3. This acts to suppresses spurious 
fragmentation since neither the Jeans mass, nor the ratio 
of Jeans length to SPH smoothing length is density depen- 
dent. Within the simulations, stars form at a rate dependent 
on their pressure. This pressure-depende nt rate is shown 
to reproduce the Kennicutt-Schmidt law (|Kennicuttlll998h 
j^t ~ A (Eg/lM0pc^)", where E* is the rate of star forma- 
tion per u nit area per unit time and is the gas surface 
density, in ISchave fc Dalla Vecchial (|20 08') . The simulations 
use a Chabrier initial mass function (jchabricr 2003) with a 
star formation rate normalisation A of 1.515 x 10~''Mq yr~^ 
kpc~^ and slope n = 1.4. Stellar particles are assigned the 
metallicity of their parent gas particle, and their subsequent 
evolution is a function of this metallicity. 
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Table 1. A list of the OWLS runs used in this analysis. We list the identifier of the simulation (as in lSchave et al.ll201(f ) and comment 
on the subgrid physics implemented in each case. Further information can be found in Section |2] 



Name 



Description 



DMONLY 

NOSNJvfOZCOOL 

REF 

WDENS 

AGN 



Dark matter only run 

No feedback, cooling assumes primordial abundances 

Weak stellar feedback with metal cooling 

Strong stellar feedback with metal cooling 

Weak stellar feedback and AGN feedback with metal cooling 




Figure 1. An example of a halo extracted from each of the simulation runs considered (top: from the 100/i~^Mpc runs at z = 0; 
bottom: from the 25/i~^Mpc runs at z = 2); an increase in brightness corresponds to an increase in surface density. In this figure the 
halo is extracted from (left to right) the dark matter only run (DMONLY); the stellar and AGN feedback (AGN) run; the strong stellar 
feedback (WDENS) run, the weak stellar feedback (REF) run and the no feedback (NOSN_NOZCOOL) run. Top: a 2 = cluster with 
M200 = S.SxlO^"* M0 and r200 = 899 h'^ kpc. Bottom: a z = 2 galaxy with M200 = 4.1x10^^ /i"^ Mq and r200 = 394 h'^ kpc. 
In each panel r2oo is shown as a white circle, while the central region 0.25r2oo is marked by a dashed black circle. It is evident that the 
baryons have a significant effect on the relaxed cluster. 



2.1.3 Feedback 



As discussed bv lWiersma et all (|2009l ). the simulations fol- 
low the timed release of mass and energy from massive stars 
[Type II supernovae (SNe) and stellar winds] and interme- 
diate mass stars (Type la SNe and asymptotic giant branch 
stars. 

In the simulations we consider, energy is injected kineti- 
c ally (stars 'kick' nearby gas part icles) using the prescription 
of lDaUa Vecchia fc Schw3 (|2008t ). Kinetic energy is injected 
locally, and the winds are not hydrodynamically decoupled. 
The efficiency of the feedback is characterised by the mass 
loading parameter, rj, and the velocity added to the nearby 
gas particles, v^. The probability for a neighbouring particle 
i to receive a 'kick' of velocity from a new stellar particle 
j is given by rirrij / ^^J^f' rrii. Typically each stellar particle 
'kicks' rj times its own mass and adds a randomly directed 
velocity, Vm, to each 'kicked' gas particle. 

The REF simulations correspond to a 'weak' feedback 
run, where = 2 and ii„ = 600kms~^ The WDENS 



run provides a more efficient form of feedback, where the 
mass loading depends on the local gas density in the fol- 
lowing way: Vw ~ 600 km s~^ (jih /0.1cm~^)^^^ and = 
2 (i;u,/600 km s~^) ^. While the same amount of the SNe 
energy is injected into the surrounding gas particles, the dis- 
tribution between mass loading and wind velocity results in 
a higher feedback efficiency. Winds in the WDENS run are 
able to remove gas from higher mass haloes more efficiently 
than the REF model. 

The final feedback run that we consider includes feed- 
back from AGN; this is, by far, the most efficient feedback 
model that we cons ider. The AGN run is i mplemented us- 
ing the method of Booth fc Schave l|2009f) which itself is 
based on that of ISpringel et al. I (|2005^ . Black holes grow 
both through merging and gas accretion (where the mini- 
mum of the Eddington and Bondi-Hoyle-Lyttleton rate is 
assumed). For star forming gas the Bondi-Hoyle rate is mul- 



tiplied by {uh/W 



to compensate for the lack 



of a cold, interstellar gas phase and the finite resolution 



Baryons and the orbital structure of haloes 5 



(|Booth fc Schavell2009l ). The black hole is assumed to grow 
as rriBH = (1 — er)iriaccr where e^, the assumed radiative ef- 
ficiency, is 0.1. Fifteen per cent of the radiated energy is 
assumed to be coupled to the surrounding medium. The 
simulation with AGN reproduces the z = observed re- 
lations between black holes a nd the mass and velocity dis- 
persion of their host galaxies (|Booth fc Schavj |2009) as well 
as the observed optic al and X-ra y properties of the groups 
in which they reside l|McCarthv e t al. 201Q ) and the steep 
drop-off in the cosmic s tar formation rate below z — 2 
jvan de Voort et al.llioill '). 

2.2 Halo sample 

Haloes within the simulation are fi rst identified usin g the 
friends-of-friends (FOF) technique (|Davis et al.lll985l ') em- 
ploying a linking length of 6 = 0.2 tim es the mean interpar- 
ticle spacing. The subfind algorithm jSpringel et al.llioOll : 
iDolag et al.l 20091 ') is then used to separate the FOF group 



into self bound structures. The main halo itself is consid- 
ered as the main subfind structure and substructures as- 
sociated with the main halo are recorded as subhaloes. The 
final halo definition considered uses a slightly modified ver- 
si on of the spherical o verdensity (SO) algorithm described 
in lLacev fc Cold (| 19941 ). A sphere is grown around the min- 
imum potential position of a halo until a specified mean 
internal density is reached (for r2oo the overdensity is 200 
times the critical density). An SO halo consists of all parti- 
cles within this sphere. 

An example of a single halo extracted from each of 
the five 100 (25) ^~^Mpc simulations at z = 0(2) is shown 
in the top (bottom) panel of Fig. [T] where an increase in 
brightness corresponds to an increase in surface mass den- 
sity. The top row shows the most massive FOF halo identi- 
fied at 2 = 0, a cluster-sized object with M200 (mass within 
r2oo) of 3.8x 10" h-^ Mq and r2oo of 899 h'^ kpc. The bot- 
tom row shows the most massive FOF halo at 2 = 2, a 
galajcy-sized object with M200 of 4.1x10^"^ h^^ Mq and r2oo 
of 394 ft-i kpc. The panels are ordered so that the efficiency 
of galaxy formation (as measured by t he central baryon frac- 
tions, see Fig.[2land fDuffv et al.llioiol ') of the simulations in- 
creases from left to right' showing the dark matter only run 
(DMONLY), the stellar and AGN feedback (AGN) run, the 
strong stellar feedback (WDENS) run, the weak stellar feed- 
back (REF) run and the no feedback (NOSN JslOZCOOL) 
run, respectively. In each panel r2oo is shown as a solid white 
circle, while the central region 0.25r2oo is marked by the 
dashed black circle. It is clear from the top panel that the 
baryons act to make the central regions of the relaxed clus- 
ter more spherical. It is interesting to note the remarkable 
similarity between the images of the cluster in the DMONLY 
and the AGN run (two top left-most images) . The galaxy at 
2 = 2 is less relaxed and the effects of the baryons are not 
obvious. 

We have selected the 50 most massive FOF haloes from 
each of the five different physics runs (discussed above) for 
this analysis. Haloes are selected from the 100 (25) h~^Mpc 
box at z = (2). At z = the mean dark matter halo mass 
of this sample is 6 x 10^'^ M©, while at z = 2 the mean 
dark matter halo mass is 7 X 10" h-^ Mq. These values 
approximately correspond to group and galaxy-scale haloes, 
respectively, so for the sake of clarity we will refer to them 



as our 'group' and 'galaxy' samples. For each halo we use 
particles associated with the main subfind group to describe 
the smooth potential of the system in order to compute the 
orbital content of the haloes. These groups contain between 
2 X 10^ and 6 x 10^ dark matter particles at z — and 
between 2 x 10* and 5 x 10^ dark matter particles at z = 
2. This corresponds to dark matter masses of between 1 x 
10" and 3 x lO" h'^ Mq at z = and between 2 x lO" 
and 4 x 10^^ Mq at z — 2. When discussing the radial 
dependence of the orbital content, we scale our results by 
the SO definition of r2oo. 

For each halo we consider its mass, dynamical state, 
spin, concentration, velocity anisotropy parameter, halo 
shape and baryon fraction. Halo properties are computed 
using the particles belonging to the main subfind group. 
We explore effect of these properties on the orbital content 
of our sample of haloes. A brief description of how these 
quantities are computed is given below. 

The dynamical state of the halo (whether it is consid- 
ered to be relaxed or not) is measured as the displacement 
of the centre of mass from the minimum potential position 
as a fraction of 7^200. If this fracti on is less than 7 per cent 
we consider the halo to be relaxed (|Neto et al.ll2007l ). Of the 
50 most massive haloes at z = (2), 27 (20) are found to be 
relaxed. 

The spin para meter of each halo is defined as in 
iBullock et all l|200ll ) as: 



A' 



J 



where J is the angular momentum within a sphere of ra- 
dius R containing mass M. The halo circular velocity V is 
defined at a radius R as, = GM/R. This spin parame - 
ter reduces to the standard spin parameter l|Peebleslll969h 
when measured at the virial radius of a truncated singular 
isothermal halo. 

We use the halo concentrations (defined as r2oo/'"s, 
where Ts is t he characteristic sc alelength of the NFW profile) 
obtained bv lDuffv et all (|2010D . 

The velocity anisotropy parameter j3 measures the pro- 
portion of radial to tangential orbits and is given by 



1 



0.54- 



(1) 



where ctj is the tangential velocity dispersion, and (t^ the 
radial velocity dispersion. A value of /3 = corresponds to 
isotropic orbits while a value of /? = 1 corresponds to purely 
radial orbits. 

In order to characterize the halo shape, we use the def - 
inition of inertia tensor given in iBaihn fc Steinmet j (|2005l ) 



The inertia tensor is diagonalised and the eigenvalues and 
eigenvectors are computed. The values of a, 6, c are defined to 
be the square roots of the eigenvalues (where a> b > c). The 
shape parameters are defined as follows: s — c/a is used as a 
measure of halo sphericity (where s — 1 for a spherical halo) 
and T = (a^ — b'^)/{a? — c?) as a measure of the triaxiality 
of the halo [where T = 1(0) for an prolate (oblate) halo] . 
Computation of the inertia tensor in a spherical region biases 
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the s hapes towards higher sph ericity; this is corrected for 
(as in lBailin fc Steinmet dboOSl ) by adopting the empirically 
motivated modified axis ratios (c/a)^^^^ = {c/a)'^^^^_^^^ and 

{b/a)^_ ^ (6/a)±._,. 

Finally, we consider the central baryon fraction of the 
halo (t he baryon to tota l mass fraction within 0.05 r2oo, 
see also iDuffv et al.ll201Gl ). The central baryon fraction ver- 
sus M200 of the OWLS haloes are shown in Fig. [5] where 
the haloes are divided into 5 mass bins, equally spaced in 
log(M2oo), error bars represent the quartile scatter. This fig- 
ure illustrates how the central baryonic mass concentration 
is affected by the strength of the different feedback mod- 
els. The left-hand figure corresponds to 2; = haloes, while 
the right-hand figure corresponds to haloes at z = 2. As 
expected, the runs with weak or no feedback have a much 
higher central baryon concentration than the stronger feed- 
back runs. The AGN run clearly has a significantly lower 
central baryonic concentration than any of the other baryon 
runs considered here. Haloes from the AGN run therefore 
appear most similar to the haloes in the DMONLY (for ex- 
ample, see the two left-most panels in Fig.[TJ. Over the range 
of halo mass probed, the central baryon fraction does not 
appear to vary significantly as a function of the halo mass 
at z = 0. At z = 2, neither of the strong feedback runs 
are mass dependent, but in the no feedback run low-mass 
haloes have higher central concentrations of baryons than 
their high-mass counterparts, while the opposite is true for 
the weak feedback run. 



3 ORBITAL CONTENT COMPUTATION 

We aim to identify the orbital content of cosmological haloes, 
to explore how the baryonic process affects the orbital con- 
tent of dark matter haloes and to link the orbital content 
of haloes with their observable and intrinsic properties. To 
this end the orbits of dark matter particles, stellar parti- 
cles and subhaloes are integrated (using a Runge-Kutta- 
Fehlberg routine) within the smooth gravitational potentials 
of the OWLS haloes. These orbits ar e then classified using 
the sp ectral classification routine of ICarpintero fc AguilaJ 
l|l998l ) into box, tube and irregular orbits allowing for a 
quantitative comparison between different simulation runs. 
Below, we summarise the key steps of this procedure. 



(as in Ijesseit et al.ll2005l : iHoffman et al.|[201Gl : [Lowing et al.l 
we have used the SCF basis functions to reconstruct 
the potential of the haloes. Density and potential are given 
as 

p(r,e,4>)= ^ k„lrr,p-nlYlm{9,(l)), (2) 

ri,Z ,m 

^{r,e,<l>)= Y. B„i™'l>„iF,„(e,0), (3) 

where n denotes the radial expansion terms and I and m 
the angular terms. There a re two commonly us ed basis 
functions: those suggested b y IClutton-Brockl (119731 ) and by 
iHernauist fc Ostriked (Il992l ). The basis set used here is con- 
structed from the latter (the code was generously provided 
by the authors) so th at the lowest or der terms represent 
the Hernquist profile (|Hernauistl Il990l ) . We find that the 
choice of basis set does not affect the reconstruction of the 
halo potentials significantly for the radial region we explore 
here (see Appendix A for further discussion). The Hernquist 
model density-potential pair is given by 

Hr)--^, (5) 
r + a 

where M is the total mass and a is a scale-length that is 
related to the half mass radius ri/2 as follows: 

a=(^-l)ri/2. (6) 

Twelve radial terms and six angular terms are used as this 
has been found to be sufficient to reproduce the potential to 
within a few percent of the A''-body potential (discussed in 
Appendix A) . 

It is important to choose a reasonable scale length in 
the Hernquist profile for the potential reconstruction. To 
optimize the potential reconstruction, particles are divided 
into two components: a diffuse component consisting of dark 
matter and hot gas (T > 10^ K) and a compact component 
consisting of stars, cold gas and black holes. The scale length 
is determined separately for each component based on its 
half-mass radius using equation ((6}, and the corresponding 
potentials for these components are computed. These poten- 
tials are then summed to give the resulting potential of the 
system as a whole. 



3.1 Calculating the potential 

The first step in determining the orbital content of a halo 
is to estimate the gravitational potential of the halo. There 
are a number of techniques that can be used to do this. As 
galaxies are regarded as coUisionless systems, an estimate 
of the smooth mean gravitational field of the system which 
minimises the effects of discrete particle representations on 
the halo potential is particularly useful. One such approach 
is the self consistent field (SCF) method. This method is 
used to obtain an estimate of the mean gravitational field 
by expanding the density and potential into a set of basis 
functions and using the forces derived from this expansion 
to integrate the equations of motion of the particles. If the 
first few terms of the basis are sufficient to provide a good 
representation of the system, then higher order terms may be 
neglected, minimizing the effects of discreteness. In this work 



3.2 Computing the orbits 

We consider the orbits of dark matter particles, stellar par- 
ticles and subhaloes. In determining the orbital content de- 
scribed by the particle distribution, a subsample of 500 par- 
ticles is selected from each halo, and the orbits of these par- 
ticles are followed within the underlying potential of the 
halo (as estimated using the SCF basis functions). One hun- 
dred particles are chosen, at random, from five radial bins 
equally spaced in log(r), with the outermost bin edges de- 
fined to be at 0.048,0.072,0.109,0.166 and 0.251 x r2oo- 
In this way we focus on the central region of the haloes 
where baryons are expected to dominate (see discussion on 
convergence testing in Appendix A). We assigned particles 
to a bin according to their 'initial' position. We have ex- 
plored binning by energy and found similar results. Par- 
ticles in the innermost region are integrated for 100 Gyr; 
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Figure 2. The baryon fraction within 5 per cent of r200 versus Af200 of the halo in each of the simulation runs. Error bars represent 
the quartile scatter. The left-hand panel corresponds to ^ = 0, and the right-hand panel to ^ = 2. The strength of the different feedback 
models is clearly visible. 



this time interval is then increased with radius. Orbits are 
also computed for all subhaloes with masses greater than 
10^°fe"^MQ. Subhaloes are integrated for 1000 Gyr. The 
motion of each particle/subhalo is integrated assuming that 
the potential remains static. A static potential is adequate 
for the purposes of this paper as we are interested only in 
characterizing the orbital content of a halo at a given point 
in time; we do not consider the evolution of this quantity. 
Also, the figure rotation of these haloes is assume d to be slow 
ijBaihn fc Steinmetj[200i : iBrvan fc Cr"es3l2007l ) and would 
probably have a negligible effect on the quantities calcu- 
lated here. A full investigation of the figure rotation of these 
haloes is deferred to future work. 



3.3 Classifying the orbit 

To classif y the orbits obtained, the s pectral classification 
routine of ICarpintero fc AguilaJ l|l998l ). hereafter CA98, is 
used. This is a fully automated classification routine based 
on the Fourier spectra of the motions of the particles. For 
a full description of the technique the reader is referred to 
their paper, but the method behind this routine is discussed 
briefiy here. 

CA98 uses the result that, once a frequency spectrum of 
an orbit is decomposed into its fundamental frequencies, the 
relationship between these frequencies can be used to classify 
the orbit in a 3D potential into the major orbital families: 
box, major(a::)-axis and minor(z)-axis tubes where orbits are 
orientated such that the major axis corresponds to the x- 
axis and the minor axis to the z-axis. Since regular orbits 
are quasi-periodic, the Fourier spectra of the time series of 
each coordinate will consist of discrete peaks (this is not 
the case for irregular orbits). The Fourier transform of the 
time series of each coordinate is performed and the dominant 
peak frequency determined. For each pair of coordinates (x- 
y, y-z and x-z) these frequencies are compared, searching 
for linear combinations (resonances). If the peak frequency 
in each direction of motion i is represented by Ui, then a 
resonance is defined as 

Iljx + mujy + nujz = 0, (7) 

for non-trivial combinations of the integers n, I and m. If 



all dominant frequencies are a multiple of a single unit fre- 
quency, then there is one base frequency. If there is no res- 
onance, then all dominant frequencies are irrationally re- 
lated. Once the dominant frequencies have been compared, 
the spectra are searched for additional base frequencies. The 
number of base frequencies specifies whether an orbit is reg- 
ular (open, closed or thin) or irregular, while the number 
of resonances specifies the orbital family as box or a;-tube, 
y-tube or z-tube. Only particles/subhaloes that have under- 
gone at least 40 orbits are classified; this ensures that they 
have clearly defined spectra. 

A 3D orbit with 4 or more base frequencies is classi- 
fied as irregular; if it has 3 (or fewer) base frequencies it is 
classified as regular. The base frequencies of a box orbit are 
incommensurable; this is the only class which does not ex- 
hibit resonance between the dominant frequencies. The orbit 
is classified as a z-tuhe if the x- and y-spectra show a 1:1 
resonance, that is I = 1, m — 1 and n is arbitrary. If y and 
z show a 1:1 resonance (m = 1 and n — 1), then the orbit 
is classified as an x-tuhe. As orbits around the intermedi- 
ate axis are unstable, it is only in rare cases that y-tubes 
are identified. These show resonances between the x- and 
z-base freq uencies. A summary of the orbit classifications 
(taken from lCarpintero fc Aguiiailll998l ) is given in Tabled 

Examples of the orbital types extracted from the dark 
matter only simulations are shown in Fig. [3] The left (right) 
panels show an example of a box (z-tube) orbit. The x-y 
projection of the orbit is shown in the top row while the 
Fourier spectra of the x- and y-motion are shown in the 
middle and bottom panels, respectively. Peaks identified by 
the routine are marked as vertical lines. Both of these or- 
bits have three base frequencies. While the box orbit has no 
resonances one can clearly see that the dominant peaks in 
the x- and y-spectra of the tube orbit show a 1:1 resonance 
{fx/ fy = 1) - this is the only resonance found. 

The CA98 algorithm has been tested rigorously using 
a number of analytic potentials. As it is fully-automated, 
it allows for the classification of large numbers of orbits. It 
also distinguishes more orbital classes than classifications 
based on the sign of a component of the orbits' angular mo- 
mentum. For comparison, the orbits considered here have 
also been cla ssified using the spin classification technique 
teamed [l993 ). While the fraction of box orbits is in gen- 
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Figure 3. Examples of the orbital types extracted from the dark matter only simulations. The left (right) columns show an example of 
a box (z-tube) orbit. The x-y projection of the orbit is shown in the top panel; the Fourier spectra of the x- and y-motion are shown in 
the middle and bottom panels, respectively. Peaks identified by the routine are marked as green vertical lines. Both of these orbits have 
three base frequencies. While the box orbit has no resonances, one can clearly see that the dominant peaks in the x- and i/-spectra of 
the tube orbit show a 1:1 resonance {fx/ fy = 1) - this is the only resonance found. Colour indicates the time evolution of the orbit, from 
red to green to blue indicates progression with time. 



eral higher than that obtained using the method of CA98, 
the same general trends are found using both classification 
schemes. 



4 RESULTS 

In this section the results of the spectral analysis of the or- 
bital content of the OWLS haloes are presented. Particular 
emphasis is placed on the fraction of box orbits, as these 
orbits are known to be important in conveying information 
from the central region to the outskirts of the halo and are 
thought to be responsible for supporting the triaxial shape 
of haloes. We quantify the fraction of different orbital types 



and show how the orbital content is affected by the addition 
of baryons and feedback processes. We then show how the 
orbits of dark matter particles are influenced by halo prop- 
erties (such as concentration, shape, spin and central baryon 
fraction) . Finally, we compare the orbits of dark matter par- 
ticles to those of stellar particles and subhaloes. 



4.1 Orbits of dark matter particles 

We begin by considering the orbital content of dark matter 
particles in each of the five simulation runs. We focus on 
the central region of the haloes (within 25 per cent of r2oo) 
where baryonic physics is likely to play a large role. These 
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results are presented as the percentage of each type of orbit 
as a function of radius in Fig. [l] For clarity we show only 
the fraction of box (black squares), loop (blue circles) and 
irregular (red triangles) orbits. The number of orbits that 
are not classified (< 10 per cent) can be determined by sub- 
tracting the sum of box, tube and irregular from 100 per 
cent. 



4-1.1 Dark matter only simulations 

The orbital content of the DMONLY haloes is shown in Fig. 
131 In the left-hand panel we show the orbital content of the 
z = haloes (mean halo mass of 6 X 10" /i"^Mq); haloes 
at a = 2 (mean halo mass of 7 x 10^^ /i~^Mq) are shown in 
the right-hand panel. Symbols show the median fraction of 
orbits of a particular type averaged over all haloes and error 
bars represent the quartile halo-to-halo scatter. In these two 
panels, we compare the complete sample of 50 haloes (solid 
lines) to the subset of these haloes that are found to be 
relaxed (dashed). The orbital content of the relaxed sample 
does not appear to differ significantly from that of the whole 
sample so all 50 haloes are used for the rest of the analysis. 

At both redshifts, the haloes are dominated by box or- 
bits out to 0.25r2oo- The dominance of box orbits is un- 
surprising; these orbits are required to support the triaxial 
haloes characteristic of dark matter simulations. The frac- 
tion of tube orbits in the galaxy sample at z = 2 is con- 
siderably lower than that found in the 2 = group sample, 
while the fraction of box orbits is slightly higher in the z — 2 
galajcy sample. We find that the fraction of irregular orbits 
in the 2 = 2 galaxy sample is much higher than that seen 
in the z — group sample, as expected when major merg- 
ers do minate the forma tion process or the mass accretion is 
rapid (|Zhao et al.ll2009l ). 

There is a weak trend for the fraction of box orbits to 
decrease with increasing radius; this is accompanied by an 
increase in the fraction of tube orbits. Resonant box orbits 
(defined in Tabled considering integer values up to n, m = 
12) account for approximately half of the box orbits shown 
here and are also found to decrease with increasing radius. 
The fraction of y-tubes, irregular and non-classified orbits 
is negligible. The fractions of both x- and z-tubes increase 
with radius. While x-tubes dominate the tube contribution 
at small radii, the fraction of 2-tubes becomes increasingly 
important at larger radii (not shown). 



4-1.2 Bar-yon simulations and the effect of feedback 

The central baryonic mass concentration is significantly af- 
fected by the strength of the different feedback models, as 
shown in Fig. [2] In this section the impact of this central con- 
centration on the orbital content of the haloes is discussed. 
The impact of baryons on the orbital content of haloes can 
be seen in Fig. [5] which is presented in the same way as Fig. 
[4] The 2 = (2) sample is shown in the left (right) column. 
In these plots the solid lines show the fraction of orbits in 
the run with baryons while the dashed lines show the orbital 
content of the dark matter only simulation for comparison. 

All of the baryon runs are found to have a smaller frac- 
tion of box orbits (at all radii out to 0.25r2oo) than the 
dark matter only haloes, but this decrease is most notice- 
able in the very central regions where baryonic condensation 
is most significant. The central concentration of baryons acts 
to transform the box orbits into tube orbits. This is a result 
of the decrease in the elongations of the orbits in response to 
the central mass (|Dubinskilll99i ). While the orbital content 
of haloes extracted from the AGN run is remarkably similar 
to the dark matter only haloes, the efficient cooling in the 
weak feedback and no feedback runs shows the most signif- 
icant reduction in the fraction of box orbits in the central 
region. These results are not unexpected. The AGN feedback 
expels most of the baryonic component from the central re- 
gions (as is evident in Fig. [2]). The runs with no or weak 
feedback have a much higher central baryon concentration 
and hence fewer box orbits than the stronger feedback runs. 
These results are in accord with expectations that increased 
galaxy formation efficiency /central baryon fraction lowers 
the fraction of box orbits. 

The right column of Fig. [5] shows the orbital content of 
the 50 most massive haloes at 2: = 2, where the mean dark 
matter halo mass is 7 x 10^^ /i~^M0. Again we note that 
there is a significant decrease (increase) in the fraction of 
tube (irregular) orbits compared to the z = Q sample. This 
is particularly apparent in the weak feedback (REF) and 
no feedback (NOSN J^OZCOOL) simulations. These runs 
appear very similar in the most central regions; perhaps un- 
surprisingly in that they appear to share very similar baryon 
fractions at z = 2 (see Fig. [2}. 

4.1.3 Orbital content versus halo properties 

In this section the dependence of the orbital content of the 
haloes extracted from the cosmological simulations on sev- 
eral key halo parameters is considered. The effect of the 
halo mass, concentration, velocity anisotropy, spin, spheric- 
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Figure 4. The orbital content of the dark matter only simulation. We show the median percentage of dark matter particles on box 
(black squares), tube (blue circles) and irregular orbits (red triangles) estimated over the 50 most massive haloes (solid lines) and for 
the relaxed subsample (dashed lines). Error bars show the quartile halo-to- halo scatter. The left-hand panel shows the orbital content of 
the haloes at z = where the mean dark matter halo mass is 6 X 10^'' h~^M0. The orbital content at z = 2 is shown in the right-hand 
panel; this sample has a mean dark matter halo mass of 7 X 10^^ ^."^Mq. 



ity and central baryon fraction on the percentage of box 
orbits can be seen in Fig. |6] As the radial dependence of the 
orbital content is found to be weak, we focus here on the 
orbital content averaged over the inner region of the haloes 
(within r < 0.25r2oo)- 

For each of the halo properties considered we divide the 
halo sample (from each run) into two subsets. The first con- 
taining the 25 haloes with the highest value of the specified 
property and the second the 25 haloes with the lowest value 
of the same property. 

Fig.[6]emphasizes the impact of baryons on the fraction 
of box orbits and illustrates their effect on several other halo 
properties. From top-left to bottom-right, the panels show 
halo mass, concentration, velocity anisotropy, spin, spheric- 
ity and central baryon fraction, respectively. In each panel 
one can clearly see that the fraction of box orbits is inversely 
proportional to the galaxy formation efficiency of the simu- 
lation. From the top-right panel it is clear that strong (weak) 
feedback runs result in haloes that are less (m ore) concen- 
tratcd than the dark matter only case (as in iDuffv et al.l 
L2OIO). It is also clear that efficient cooling results in more 
spherical haloes (bottom-left panel). 

From Fig.[6]we can see that the percentage of box orbits 
is not sensitive to the halo mass, concentration, anisotropy 
or spin parameter (although there is a tentative trend for 
haloes with high spin parameters to have fewer box orbits 
in the weak/no feedback runs) over the range of parame- 
ters considered here. Trends are apparent when we consider 
halo shape and central baryon fraction. It is clear from the 
bottom-left panel of Fig. |6] that an increase in sphericity 
corresponds to a decrease in the percentage of box orbits. 
It is also clear from the bottom-right panel that the cen- 
tral baryon fraction has a significant effect on the orbital 
content; an increase in central baryon fraction corresponds 
directly to a decrease in the fraction of box orbits. 

4.2 Orbits of stellar particles 

Next we compare the orbits of dark matter particles to those 
of the stellar particles. We note, however, that since the po- 
tential and initial conditions are drawn from the main sub- 



find halo, distinct subhaloes and satellites are not included 
and the stellar population considered here is associated with 
the central galaxy and the diffuse intra-halo component. 

We show the orbital classifications of stellar particles in 
Fig. [T] presented in the same way as in Fig. U) Symbols show 
the median fraction of orbits of a particular type averaged 
over all haloes and error bars represent the quartile halo- 
to-halo scatter. The dashed curves show the orbits of the 
dark matter particles taken from the same simulation, for 
comparison. In the left (right) column we show the orbital 
content of the 2 = (2) haloes. As in Fig. 31 the different 
feedback implementations are compared, with the galaxy 
formation efficiency of the simulations increasing from top 
to bottom. 

The orbital content described by the stellar particles at 
2; = is remarkably similar to that drawn from the orbits 
of dark matter particles. While the dark matter and stellar 
particles are selected from the same radius, one might ex- 
pect a different trend due to the formation history of the 
stellar particles. A full analysis of the history of the stellar 
particles (such as when they were stripped from parent sub- 
haloes) could prove insightful. The results presented here 
seem to indicate that they were either stripped a long time 
ago and have forgotten their dynamical history, or subhaloes 
bringing in stars are not biased significantly with respect 
to the main d istribution - the velocity bias is weak (e.g. 
ISpringel et "ali r2001 ). 

At 2 = 2 the same trends seen in the orbital content of 
the dark matter particles are visible in the orbital content 
of the stellar particles. However, the fraction of stellar parti- 
cles on irregular orbits is significantly enhanced compared to 
that for dark matter particles, coming at the expense of the 
box orbits. The fraction of dark matter and stellar particles 
on tube orbits is remarkably similar. 

As a final comparison between the orbits of dark matter 
and stellar particles we show, in Fig. [8] the dependence of the 
fraction of stellar particles on box orbits (within 0.25r2oo) 
on the halo mass and central baryon fraction (within 0.05 
7'20o) • This is directly comparable to the top-left and bottom- 
right panels of Fig. [6] Once again we see a clear indication 
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Figure 5. The median percentage of dark matter particles on box (blaclc squares), tube (blue circles) and irregular orbits (red triangles) 
estimated over the 50 most massive haloes in each simulation (solid lines). Error bars show the quartile halo-to-halo scatter. Dashed lines 
show the orbital content of the 50 haloes in the dark matter only run, for comparison. The left-hand panels show the orbital content of 
the haloes at z = 0, where the mean dark matter halo mass is 6 X 10^^ /i^^Mq. The orbital content at 2 = 2 is shown in the right-hand 
panels (mean dark matter halo mass of 7 X 10^^ h~^'MQ). 
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Figure 6. The effect of basic iialo properties on the fraction of box orbits in the central region (within 0.25r2oo) of group-sized haloes 
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the specified property and the second the 25 haloes with the lowest value of the same property. Error bars show the quartile halo-to-halo 
scatter for the subset of haloes. 



of the impact that baryons have on reducing the fraction of 
box orbits. 



4.3 Orbits of subhaloes 

Finally, we consider the orbits of the subhaloes associated 
with each SUBFIND main halo at 2 = 0. We consider all 
subhaloes with masses greater than 10^" /i~^M0, there are 
not enough objects to restrict the sample to the central (r < 
0.25r2oo) region. The subhaloes chosen in this way trace a 
region much further out than discussed previously (mean 
r/r2oo = 0.6 compared to a mean value of r/r2oo = 0.12 for 
the dark matter and stellar particles) and are likely to probe 



a region less strongly affected by the presence of baryons. 
Subhaloes on box orbits are also likely to be strongly affected 
by tidal disruption. 

The initial position of each subhalo was taken to be an 
average over the 10 most bound particles, and the velocity of 
the subhalo is assumed to be that of the most bound particle. 
The orbits of these subhaloes were integrated, within the 
gravitational potential of the main halo, for 1000 Gyr. After 
this ~85 — 90 per cent of subhaloes from the baryon runs 
are classified; however, only ~70 per cent of the dark matter 
subhaloes had undergone more than 40 orbits. 

The fraction of subhaloes on box orbits as a function 
of the halo mass M200 and central baryon fraction (within 
0.05 r2oo) is shown in Fig. (9] All of the baryon runs indicate 
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Figure 7. The median percentage of stellar particles on box (black squares), tube (blue circles) and irregular orbits (red triangles) 
estimated over the 50 most massive haloes in each simulation (solid lines). Error bars show the quartile halo-to-halo scatter. For 
comparison, in all panels, dashed lines show the percentage of dark matter particles (from the same simulation) on each type of orbit. 
The left-hand panels show the orbital content of the haloes at ^ = 0, where the mean dark matter halo mass is 6 X 10^^ /i^^Mq. The 
orbital content at 2 = 2 is shown in the right-hand panels (mean dark matter halo mass of 7 X 10^^ h~^MQ). 
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Figure 8. The fraction of stellar particles on box orbits for the 
different simulation runs at z = 0. Top: fraction of box orbits as a 
function of the halo mass A'/2oo- Bottom: fraction of box orbits as 
a function of the central baryon fraction (within 0.05 r2oo)- Error 
bars show the quartile halo-to-halo scatter. 



Figure 9. The fraction of subhaloes on box orbits for the dif- 
ferent simulation runs at 2 = 0. Top: fraction of box orbits as a 
function of the halo mass M200 • Bottom: fraction of box orbits as 
a function of the central baryon fraction (within 0.05 r2oo)- Error 
bars show the quartile halo-to-halo scatter. 



a similar fraction of box orbits; this fraction is higher than 
that seen in the dark matter only simulation only because 
of the fraction of subhaloes that remain unclassified in this 
simulation. If we consider the fraction of classified subhaloes 
on box orbits, then the fraction in the dark matter only 
simulation increases to ~60 per cent, in agreement with that 
seen in the AGN run. 

Comparing the orbital content of the subhaloes to that 
of the dark matter particles (Fig. |4] and Fig. [5}, we find that 
the orbits of the subhaloes are in broad agreement with those 
seen in the outermost radial bins of the particle distribu- 
tions. Perhaps suggesting that subhaloes bringing stars into 
the main galaxy are not biased significantly with respect to 
the main distribution and explaining why the orbits of the 
diffuse intra-halo stellar component are so similar to those 
of the dark matter particles. 



5 SUMMARY AND DISCUSSION 

The orbital content of a large sample of haloes extracted 
from state-of-the-art high-resolution cosmological hydrody- 
namical simulations from the OWLS project is explored in 
order to identify potential signatures of the formation pro- 
cess. We focus on the central regions (0.25r2oo) of haloes 
with virial masses ~ 6 x 10"(~ 7 x lO") H'^Mq at 2 = 
(2) and study how the orbital content of these haloes is af- 
fected in the presence of baryons. 



Haloes in dark matter only simulations are dominated 
by box orbits out to 0.25r2oo- This is not surprising as box 
orbits are known to be required to support the triaxial haloes 
characteristic of dark matter only haloes and are thought to 
dominate in systems that have undergone major mergers. 
The fraction of box orbits is found to decrease with increas- 
ing distance from the halo centre; this is mirrored by an 
increase in the fraction of tube orbits. While at 2 = very 
few of the orbits were classified as irregular, at 2 = 2 irreg- 
ular orbits are more common. 

While spherical haloes tend to have fewer box orbits, the 
orbital content of the haloes does not appear to be strongly 
dependent on halo properties such as mass, concentration, 
velocity anisotropy, spin and dynamical state for the range 
of parameters considered here. It is, however, strongly de- 
pendent on the central baryon fraction. 

By comparing simulations run with no feedback, with 
stellar feedback and with feedback from AGN, the fraction 
of box orbits in the central region is found to decrease when 
baryonic physics is included. Baryons are able to cool and 
condense to the centre of the halo, and this central con- 
centration tends to transform box orbits into tube orbits. 
Increasing the strength of the feedback implementation is 
found to reduce the central concentration of baryons and 
increase the fraction of box orbits. The orbital content of 
the strongest feedback run (AGN) is very similar to that 
seen in the dark matter only case. 

We then compared the orbital content described by the 



Baryons and the orbital structure of haloes 15 



dark matter particles to that of the stellar distribution and 
the subhaloes. The orbital content described by the stellar 
particles (within the central galaxy and the diffuse intra-halo 
light) is found to be remarkably similar to that drawn from 
the orbits of dark matter particles. Typically ~ 50 - 60 per 
cent of the subhaloes are found to be on box orbits regard- 
less of the baryonic physics implemented. This fraction is 
in broad agreement with that found in the outermost radial 
orbits of the dark matter particles. The subhaloes we have 
analysed probe a more extended region of the halo where 
the effects of baryons do not appear to be as significant. 

While the results presented here highlight the impor- 
tance of the baryons on the orbital content of haloes, we 
are limited by the resolution of the simulations. Ideally this 
analysis would be extended to a study of the innermost re- 
gions of galaxies. The stellar half-mass radius of the haloes 
considered here is well below the innermost radial bin that 
we are able to consider. By studying high-resolution haloes 
that have been resimulated from cosmological conditions, 
we would be better placed to make direct comparisons with 
observations. The Gaia satellite, soon to be launched, will 
provide us with a kinematic census of our Galaxy and place 
strong constraints on galaxy formation models. In order to 
fully exploit such observational datasets, a comprehensive 
comparison with simulations is essential. This has moti- 
vated much work o n the topic including a recent paper by 
IValluri et al] l|201lD who find, in agreement with our study, 
that orbital analysis can provide constraints on the under- 
lying potential and may prove to be a useful probe of the 
formation history of the system. 

Analysis of the sort presented here may also prove 
useful to modelling approache s such as Schwarzschild's 
method (ISchwarzschildl 1979^, Made-to-Measure tech- 



niques (ISver fc Tr emaind Il996l: Ide Lorenzi et al. 1 l2007l : 
ao| 



iDehnenl l2009l: [Lon g & Ma^ 1201 d ) and torus modelling 



I McMillan fc BinnevI 120081 ). These methods attempt to re- 
produce observables using a superposition of the orbital 
density distributions. While these techniques are, in many 
cases, able to reproduce a number of observables, the cho- 
sen orbital distribution may not be unique. It is possible 
that different combinations of orbits with distinctly differ- 
ent s hapes may produce the s ame triaxial density distribu- 
tion l|de Zeeuw fc Fraru3ll99ll ). Studying the orbital content 
of simulated merger remnants provides unique insight into 
the types of orbital distributions expected to be present in 
galaxies with specific properties, and may help to provide 
initial conditions, or additional constraints, for these meth- 
ods. 

With better future data (particularly those from inte- 
gral field units), a comparison of orbit classifications from 
dynamical models and cosmological simulations may provide 
insights into baryonic processes in galaxy formation. 
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APPENDIX A: NUMERICAL ISSUES 

In this section the convergence radius rconv and resolution 
effects are considered in order to show that the results pre- 
sented in this paper are well converged. The choices of halo 
definition and basis sets are also discussed. 

Al Convergence radius 

iPower et al.] (|2003l ) showed that numerical convergence in 
the inner regions of dark matter haloes was achieved out- 
side of the convergence radius, Vconv, defined to ensure that 
the two-body dynamical relaxation time within this radius 
is comparable to the age of the Universe. The convergence 
radius depends on halo size and on the resolution of the 
simulation; it sets a minimum resolved length scale for the 
analysis. Fig. lAll shows an estimate of the convergence ra- 
dius as a fraction of M200 for the haloes used in this analysis. 
The horizontal line shows the innermost radial bin we con- 
sider; this has been chosen to ensure that our results are 
converged. 

A2 Resolution effects 

In order to quantify the effects of resolution on the orbital 
content, the 512^-particle run from the DMONLY simula- 
tions (with a maximum softening length of 2 kpc) is 
compared with the corresponding lower resolution runs (con- 
taining 256'^ and 128'^ particles and with maximum soften- 
ing lengths of 4 and 8 kpc, respectively). In Fig. IA2l the 
fraction of box orbits found in the 128'^ simulation is shown 
in blue, the 256^^ simulation in green and the 512^ simula- 
tion used for this analysis in red. The fraction of box orbits 
found in the dark matter only run is shown in the top plot 
while the fraction of box orbits found in the weak feedback 
run (REP) is shown on the bottom. Only relaxed haloes 
that are matched between the different resolution runs are 
considered. Vertical lines show the convergence radii for each 
simulation; bins that are considered to be converged are con- 
nected by solid lines to aid in the comparison. 

A3 Effect of halo definition 

To explore the effect of the halo definition on the results 
presented here, the three common definitions of a group - 



Baryons and the orbital structure of haloes 17 



0.045 
0.04 



^ 0.03 
J 0.025 



0.015 



0.01 



▼ 

▼ 


DMONLY T 








'V 








■» » T 








w 

T 



100 



2.8 



13 



13.2 13.4 13.6 
log M^oo ( h'^ MJ 



13.8 



14 



14.2 



0.045 p 

0.04 
0.035 

r 

0.03 
J 0.025 ' 
J 0.02 
0.015 
0.01 
0.005 

<)o^ 



11 11.2 11.4 11.6 11.8 12 12.2 12.4 12.6 
1 I 



Figure Al. Convergence radius Tconv for the dark matter only 
simulation at 2 = (top) and 2 = 2 (bottom). The horizontal 
line shows the innermost bin considered in this analysis. Orbits 
of particles beyond the convergence radius are studied. 



FOF, main subhalo as identified by SUBFIND and the SO ap- 
proach - were used. The orbital content of the haloes is not 
significantly affected by the halo definition, as shown in Fig. 
IA3I In the dark matter only simulations the fraction of box 
orbits does not depend on the choice of groupfinder. In the 
weak feedback run (REF) the subfind main haloes show a 
slightly lower fraction of box orbits due to the effect that 
removing subhaloes has on the halo potential. This differ- 
ence is not sufficient to account for the trends discussed in 
this paper. The main subhalo is therefore used throughout 
this analysis. This has the advantage of providing a smooth 
potential, unperturbed by substructure. 
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Figure A2. The effect of resolution on z = results. Results 
from the 128"^ simulation are shown in blue, the 256^ simulation 
in green and the 512'^ simulation used for this analysis in red. 
The fraction of box orbits found in the dark matter only run is 
shown in the top panel. The fraction of box orbits for dark matter 
particles found in the weak stellar feedback run (REF) is shown 
in the bottom panel. Error bars show the quartile halo-to-halo 
scatter. The convergence radius for each simulation is shown as a 
vertical line. 



bital classifications at the level of a few percent; a similar 
orbital content is found using both basis functions. 

The choice of basis set is not found to affect the recon- 
struction of the potential significantly over the radial range 
we consider here. Both basis sets reproduce the A''-body po- 
tential with percent level accuracy over this region (see Fig. 

a. 



A4 Choice of basis set and expansion coefflcients 

The basis set used in this analysis is constructed so that 
the lowest orde r terms represent the Hernquist profile 
l|Hernauistl 19901 ). Twelve radial terms and six angular terms 
are used as this has been found to be sufficient to reproduce 
the potential to within a few percent of the A''-body poten- 
tial (illustrated in Fig |A4p . 

Table lAll shows the orbital content of 500 dark mat- 
ter particles taken from the inner region (0.25r2oo) of the 
most massive relaxed cluster from the weak feedback sim- 
ulation at z = 0. The orbital content as determined using 
the Hernquist basis set for a number of different expansion 
coefficients is shown in the top panel. The orbital content 
derived using the Glutton-Brock basis sets is shown below. 
Varying the number of expansion coefficients affects the or- 



18 S.E. Bryan et al. 



Table Al. Classifications of the orbits of 500 dark matter particles taken from the inner region (0.25r20o) of the most massive relaxed 
cluster from the weak feedback simulation at 2; = 0. Different numbers of expansion coefficients have been used to reconstruct the 
potential using the Hernquist basis set (top) and the Clutton-Brock basis set (bottom). 
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Figure A4. Difference in potential as computed by the SCF ba- 
sis functions (using n = 12, I = 6) and the Direct-Summation 
approach for the most massive weak feedback (REF) halo. For 
the green points the potential has been computed using the basis 
set of Hernquist-Ostriker, while for blue triangles the potential is 
calculated using the Clutton-Brock basis set. 
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Figure A3. The effect of halo definition on the orbital content 
of haloes, using the 512^-particle simulations at ^ = 0. The top 
plot shows the fraction of box orbits in the dark matter only 
simulations, while the bottom plot shows the fraction of box orbits 
in the weak feedback run (REF). Results obtained if haloes are 
defined using the FOF algorithm are shown in red, those using the 
main SUBFIND halo are shown in green and those obtained using 
SO are shown in blue. Error bars show the quartile halo-to-halo 
scatter. The orbital content of a halo is not particularly sensitive 
to the halo definition. 



